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Introduction 

The instabilities responsible for the first two 
stages of transition to turbulence in wall-bounded 
shear flows are now quantitatively well understood. 
The first, or linear, stage is dominated by slowly 
growing Tollmien-Schlichting waves propagating in 
the streamwise direction, and the development of 
these waves is described quite accurately by the clas- 
sical linear theory. The landmark experiment of 
Schubauer and Skramstad (1947) provided the con- 
firmation of Tollmien-Schlichting waves in the Bla- 
sius boundary layer. Similar demonstrations for con- 
trolled boundary layers and for channel flow have 
been furnished by Strazisar, Reshotko, and Prahl 
(1977) and by Nishioka, Asai, and Iida (1980), 
respectively. 

The second stage is characterized by the in- 
stability of finite amplitude, streamwise (or two- 
dimensional) Tollmien-Schlichting waves to infinites- 
imal, oblique (or three-dimensional) disturbances. 
Orszag and Patera (1983) and Herbert (1983a, 
1983b) have used linear Floquet theory (see, for ex- 
ample, Coddington and Levinson 1955) to model 
these secondary instabilities. These instabilities in- 
duce strong vortical patterns in the boundary layer, 
the basic element of which is a lambda vortex. A 
regular pattern develops. It may be either ordered or 
staggered (see fig. 1). In the former case, the tips of 
successive lambda vortices are aligned in the stream- 
wise direction of flow, and the distance between 
the tips corresponds to the streamwise wavelength 
of the finite amplitude, two-dimensional Tollmien- 
Schlichting (TS) wave. In the latter case, the tips 
of successive lambda vortices are offset in the span- 
wise direction by half the characteristic spanwise 
scale of the secondary instability. Thus, the stream- 
wise wavelength of the staggered array of vortices 
is double the wavelength of the two-dimensional TS 
wave. The ordered pattern is a consequence of a 
fundamental instability, whereas the staggered pat- 
tern results from a subharmonic instability. The 
fundamental instability is called the K-type insta- 
bility after Klebanoff, Tidstrom, and Sargent (1962). 
Craik (1971) and Herbert (1983b) have identified two 
classes of subharmonic instability, known as the C- 
type and the H-type, respectively. All three classes 
have been clearly observed in experiments (Saric, 
Kozlov, and Levchenko 1984). Floquet theory pro- 
duces accurate quantitative predictions of the early 
(linear) development of these three-dimensional sec- 
ondary instabilities. 

Subsequent stages in the transition to turbulence 
are as yet understood only qualitatively if at all. The 
lambda vortices rise away from the wall and intensify 


rapidly into hairpin vortices with attendant highly 
localized, strong, detached shear layers. These regis- 
ter experimentally as one, then two, then three, etc., 
spikes in the velocity signals. After roughly the five- 
spike stage, turbulent spots occur and then fully de- 
veloped turbulence ensues. Although the instability 
of a two-dimensional, unsteady, detached shear layer 
has been modeled qualitatively by Greenspan and 
Benney (1963), the actual shear layers are decidedly 
three-dimensional. The detailed flow field structure 
at these stages is difficult to ascertain experimen- 
tally because of the very small spatial scales involved, 
not only in the turbulent region but also near the 
sharp flow field gradients connected with the hairpin 
vortices. 

Three-dimensional, nonlinear numerical calcula- 
tions of the time-dependent Navier-Stokes equations 
have simulated many important features of transition 
and turbulence in wall-bounded shear flows. Such 
simulations have thus far resorted to the parallel 
flow assumption, in which the spatial growth of the 
boundary layer is ignored, in order to relax the ex- 
treme resolution demands in the streamwise direction 
that arise for a true, spatially developing flow. Both 
qualitative and quantitative agreement has been ob- 
tained between these full Navier-Stokes calculations 
and the Floquet theory predictions for the secondary 
instabilities of the K-type (Orszag and Patera 1983; 
Kleiser and Schumann 1984; and Spalart 1985) and 
of the H-Type (Spalart 1985; Zang and Hussaini 
1985a; and Laurien 1986). Moreover, Wray and 
Hussaini (1984) demonstrated that numerical simu- 
lations based on the parallel flow assumption achieve 
good agreement with the detailed flow field structure 
of the Kovasznay, Komoda, and Vasudeva (1962) K- 
type boundary-layer experiment up to the two-spike 
stage. Likewise, Kleiser and Schumann (1984) repli- 
cated many salient features of the channel flow ex- 
periment of Nishioka, Asai, and Iida (1980). 

Even with the parallel flow assumption, neither 
of these calculations possessed sufficient resolution 
to track reliably the details of the flow all the way 
through transition. Both calculations were ceased 
once the computed solution evinced appreciable os- 
cillations on the scale of the computational grid. The 
danger, of course, in pushing a calculation beyond its 
resolution limits is that numerical artifacts may be 
mistaken for real physical effects. Orszag and Kells 
(1980) have characterized the artificial oscillations 
that arise in poorly resolved flows as “numerical tur- 
bulence.” In problems such as Rayleigh- Benard flow, 
for which bifurcation theory is reasonably successful, 
numerical turbulence may lead one to conclude er- 
roneously that a given flow is chaotic (Curry et al. 
1984). In transition problems, it could lead one to 


mislocate the instant and region in which turbulence 
originates. 

This paper has two objectives. The first one is 
to establish firmly, both for this and for subsequent 
work, the resolution requirements at various stages 
in the transition process. A single transition simu- 
lation is examined in detail. Comparisons are made 
between second-order finite-difference methods and 
spectral methods and between calculations made on 
different grids. The latter comparisons suggest a 
set of guidelines for when and in which directions 
grid refinements are required to maintain a physi- 
cally meaningful calculation. 

The second objective is to map out, at the latest 
stage currently feasible, the structure of the hairpin 
vortices that arise in K-type and H-type transitions 
in channel flow. Zang and Hussaini (1985b) have dis- 
played this vortex structure at the one-spike stage for 
parallel boundary-layer transitions, including a case 
with wall heating control. Their results agreed with 
the general features measured by Williams, Fasel and 
Hama (1984). The focus here is on channel flow, in 
part to furnish a comparison between the channel and 
boundary-layer structures, but principally because 
current spectral numerical algorithms permit more 
efficient and accurate simulation of channel flow. 

This investigation has been conducted for channel 
flow at a subcritical Reynolds number and at a su- 
percritical Reynolds number. Secondary instabilities 
of both the K-type (Orszag and Patera 1983) and the 
H-type (Zang and Hussaini 1985a) have been exhib- 
ited under the subcritical conditions. 

Although the present work is focused on the hair- 
pin vortices that arise in the transition process, it 
is also relevant to the vortical structures that have 
been observed in turbulent boundary layers. Wal- 
lace (1982) has noted that the properties of the hair- 
pin vortices in transitional and turbulent boundary- 
layer flows are similar. Smith (1984) has measured in 
detail hairpin vortices in turbulent boundary layers. 
Furthermore, several numerical simulations of tur- 
bulent channel flow have mapped out hairpin vortex 
structures (Moin and Kim 1982; and Moser, Moin, 
and Leonard 1983). 

Formulation 

The channel flow geometry is illustrated in fig- 
ure 2. Lengths are scaled by the channel half-width 
h, velocities by the centerline velocity u c , and the 
density by its constant value. The Reynolds num- 
ber R is based on these scales; that is, R = u c hfv , 
where v is the kinematic viscosity. The flow is as- 
sumed to be fullv developed with the parabolic profile 
uo(F) = 1 — Y . The incompressible Navier-Stokes 


equations for the dimensionless variables are used in 
the form 

u^o;xu = - VP -f- -)-V 2 u (- 1 ) 

R 

V • u = 0 (2) 

with no-slip boundary conditions at the walls 
(y = ±1) and periodic boundary conditions in X 
and Z. The velocity vector is denoted by u, the vor- 
ticity a? = Vxu, and the total pressure P = p-b ^ |u| 2 , 
where p is the static pressure. 

Two related numerical methods have been used to 
obtain the results reported herein. Both use a Fourier 
collocation technique for the spatial discretization 
in the X- and Z-directions. One uses a second- 
order finite-difference scheme in Y, and the other 
uses Chebyshev collocation for that direction. The 
Fourier finite-difference method (FFDM) was devised 
by Moin and Kim (1982) for their large-eddy simula- 
tion and was used by Biringen (1984) for his transi- 
tion calculations. The particular Fourier-Chebyshev 
method (FCBM) used here is the one developed by 
Zang and Hussaini (1986). The FCBM code uses the 
standard Chebyshev grid in Y, whereas the FFDM 
code uses the hyperbolic tangent transformation, 

y = itanh(^tanh _1 a) (3) 

The computational variable £ € [ — 1 , 1] , and the 
stretching parameter a = 0.94175 is used. Both 
methods employ an explicit time discretization for 
the advection and horizontal diffusion terms and an 
implicit time discretization for the vertical diffusion 
and pressure gradient terms as well as the continuity 
equation. The FCBM uses a splitting, or fractional 
step, technique for the time advancement: in the first 
half-step the advection and diffusion terms are ac- 
counted for, and in the second half-step a pressure 
correction is applied so that the final velocity field 
is divergence free. Boundary conditions are applied 
at the intermediate step to all three velocity compo- 
nents and for the pressure correction to the normal 
velocity component. The FFDM is unsplit; that is, 
all terms are accounted for simultaneously in a sin- 
gle step. Both codes employ a staggered grid for the 
pressure so that no artificial boundary conditions are 
required. 

There are at least five other fully spectral algo- 
rithms for channel flow simulation — see Orszag and 
Kells (1980), Kleiser and Schumann (1984), Moin 
and Kim (1980), Moser, Moin, and Leonard (1983), 
and Zang and Hussaini (1985a). Only the first 
of these resorts to splitting. Zang and Hussaini’s 
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(1985a) method is the unsplit counterpart to the 
FCBM used in the present work. The split FCBM 
code is more robust and requires less storage and 
CPU time than the unsplit code. Zang and Hussaini 
(1986) have confirmed that when appropriate inter- 
mediate boundary conditions are employed, the split 
FCBM code gives the same results (to at least four 
significant digits) as the unsplit version under the 
precise physical conditions simulated in the present 
work. 

In these algorithms, both physical space and spec- 
tral space representations of the dependent variables 
are employed. The velocity has the Fourier series 
representation, 


N x / 2-1 N g / 2-1 

<x,Y,z,t)= E 

k x = -N x / 2 k z = -N z /2 
x g* [(fc*/s*)otX+ (kz/sz}pZ^ 

where a and /? are the fundamental wave numbers in 
the streamwise and spanwise directions, respectively, 
and s x and s z denote the number of subharmonic 
modes which are included in these directions. The 
fundamental wavelengths in the X- and Z-directions 
are L x = 2 ir/a and L z — 27t// 3, and the computa- 
tional domains have lengths s x L x and s z L z . The in- 
dividual Fourier modes are labeled by (k X) k z ), where 
k x = k x /s x and k z = k z /s z . Several of these Fourier 
components are illustrated in figure 3. 

The initial conditions for the simulations re- 
ported in this work consist of the mean flow plus 
a two-dimensional TS wave and two oblique three- 
dimensional waves: 

u(X,Y,Z,0) = Re{u 0 (y) + £ 2 DU2DCne lQX 
+ ^3DuJ D (y)e* l( “ /s3:)X+,3Z! 

+ \iZDVz D {Y)e i[(a/Sx)X ~ l3Z] } (5) 

where U2£>(T) and u±p(Y) are the least stable lin- 
ear modes for the given real wave numbers a and 
f3. These eigenfunctions are normalized so that their 
maximum streamwise amplitudes are 1. Thus, 62 d 
and 63 d measure the amplitudes of the initial stream- 
wise velocity perturbations. For the K-type transi- 
tion simulation the three-dimensional waves are TS 
waves and s x = 1, whereas for the H-type transition 
the three-dimensional waves are Squire modes, that 
is, eigensolutions of the homogeneous vertical vortic- 
ity equation, and s x = 2. 


A useful measure of the strength of a given Fourier 
harmonic is the quantity 


E kx , k ,(t) = i|, d kx d kz \u kx , k ,(Y, t ) | 2 dY (6) 


where 

This quantity is the kinetic energy of the mode 
normalized by the kinetic energy of the mean flow. 

Similarly, the velocity has the Chebyshev series 
representation, 


Ny 

u(X, y, Z,t)=Y 2 Z, t)T n (Y) (7) 

n— 0 


where T n (Y) is the Chebyshev polynomial of degree 
n. The significance of these modes is measured by 


E n {t) = 


15 


16 Sx Lxs z L z 



|u n (X, Z,t)\ 2 dXdZ 


( 8 ) 

The calculations were performed in 32-bit arith- 
metic on the Control Data Corporation VPS 32 at 
NASA Langley Research Center. The FFDM calcu- 
lations took 0.65 sec per step on a 32 x 50 x 32 grid 
and 2.5 sec per step on a 32 x 200 x 32 grid. 
The FCBM code requires 0.96 sec per step on a 
32 x 64 x 32 grid and 2.7 sec, 7.50 sec, and 15.5 sec 
per step on 64 3 , 96 3 , and 128 3 grids, respectively. 
The largest grid used in the calculations reported 
here was 96 x 128 x 216. This consumed 20.1 sec per 
time step. 


Resolution Study 

The particular subcritical K-type transition stud- 
ied by Biringen (1984) is used here to illustrate the 
resolution demands of transition simulations. The 
Reynolds number is 1500. The horizontal wave 
numbers a and (3 are both unity. No subharmon- 
ics are included (s x = s 2 = 1). The temporal 
frequencies of the least rapidly damped modes are 
UJ 2D = 0.326299 - 0.028206i and c u 3D = 0.401293 - 
0.028230i, as determined by a Chebyshev tau approx- 
imation to the Orr-Sommerfeld equation (Orszag 
1971b). The initial amplitudes of the waves are 
62 d ~ 0.11 and €30 = 0.10. The primary two- 
dimensional wave has the period T = 19.26. The 
numerical computations are made in a frame with a 
streamwise velocity of 0.36. This is the average of 
the phase speeds of the two- and three-dimensional 
waves. The use of a moving frame improves the tem- 
poral accuracy of the calculation. 
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A straightforward test of the vertical resolution 
is the ability of a code to reproduce the time evolu- 
tion of a linear eigenfunction. The computed growth 
rates and their errors at t = 30 are given in table I for 
several calculations. The initial condition for each 
case consisted of a three-dimensional TS wave at 
0.01-percent amplitude. The Chebyshev calcula- 
tion is accurate to better than five significant dig- 
its. The finite-difference calculations exhibit the ex- 
pected second-order accuracy. Moreover, the finite- 
difference calculations are on the side of excessive 
damping of the wave. (This same effect was found 
by Malik, Zang, and Hussaini (1985) for a higher 
Reynolds number.) 

The error in the 32 x 50 x 32 FFDM growth rate 
might appear to be small. Nevertheless, it indicates 
uncomfortably large errors in the full transition sim- 
ulation. These manifest themselves quite clearly in 
the time-dependent vortical flow field. Since the vor- 
ticity is obtained from the primitive variables by dif- 
ferentiation, it is a sensitive indicator of errors in 
the velocities themselves. Typical vorticity patterns 
are given in figure 4, 1 * which shows the full three- 
dimensional fields in the lower half of the channel. 
(The features in the upper half are similar, but are 
displaced in X by L x f 2.) These color plots show 
the streamwise = dw/dY — dv/dZ) and span- 
wise (uj z = dv/dX - du/dY ) vorticity components 
in terms of two level surfaces each. To improve per- 
ception, the surfaces have been shaded to mimic the 
effect of light sources. Most of the hairpin vortex is 
traced out by the uppermost red-blue pair of surfaces 
in the u x plot. The tip, of course, is missing because 
c o x vanishes there. The peak plane is the center of the 
spanwise (Z) direction. It lies between the two legs 
of the hairpin vortex. The uppermost surface in the 
t o z plot is in the region of the detached shear layer. 
Note that in the spanwise direction, this shear layer 
is confined to the vicinity of the peak plane. For this 
resolution study, we focus on contour plots of a few 
key cross sections, including the peak plane (Z = 7r). 

Errors in the 32 x 50 x 32 FFDM calculation are 
reflected by oscillations in the contours of vertical 
shear (du/dY) in the peak plane even at the early 
time t = 7.5. These contours are plotted in figure 5, 
and the results from a 32 x 200 x 32 calculation are 
given for comparison. The coarse grid results are al- 
ready exhibiting the grid-scale oscillations which, in 
what is still laminar flow, are indicative of inadequate 
resolution. Figures 6 and 7 furnish a similar compar- 
ison at t = 15. The coarse grid result not only fails to 


1 The authors are grateful to Robert Weston, Langley 

Research Center, for supplying the three-dimensional color 

figures. 


pick up the full intensity of the detached shear layer 
but also misrepresents its speed. The fine grid result 
places the center of the shear layer at X & 0.55 L x , 
whereas it has reached only X as 0.45L X in the coarse 
grid calculation. 

These figures also include a (Z-Y) contour plot 
of the streamwise vorticity (u; x ) at a streamwise 
station near the maximum of this vorticity compo- 
nent (roughly 1/2 of the distance along the x-axis 
at t = 15). (The features corresponding to the 
lambda-hairpin vortex are those whose centers are 
near Y = —0.625.) The vertical shear plots are in- 
dicative of the resolution in X and F, whereas the 
streamwise vorticity plots highlight the Y and Z reso- 
lution. Even for the fine vertical grid FFDM calcula- 
tions, there are noticeable oscillations in the stream- 
wise vorticity. These are removed with increased 
spanwise resolution, as is demonstrated below. 

The fully spectral FCBM result on a 32 x 64 x 32 
grid at t = 15 is given in the top half of figure 8. 
It agrees well with the 32 x 200 x 32 FFDM result. 
The vertical shear plots differ only near the center of 
the detached shear layer, and the FFDM streamwise 
vorticity plots are slightly more oscillatory. Thus, the 
poor results of the 32 x 50 x 32 FFDM calculation at 
t = 7.5 and t = 15 are due primarily to inadequate 
resolution in the vertical direction and not to any 
flaw in the code itself. 

The results in the bottom half of figure 8 are for 
t = 18.75. The spanwise fluctuations in the stream- 
wise vorticity have the unmistakable signature of a 
grid-scale oscillation with wavelength 2 A Z. They 
evidently arise because the number of points in this 
direction is insufficient to resolve the sharp gradients 
in the vortex loop. 

As is illustrated in figure 9, a 64 3 FCBM calcula- 
tion eliminates the oscillations at t — 15, but is still 
too coarse at t — 18.75. Note the shorter scale (com- 
pared with fig. 8) of the fluctuations in the stream- 
wise vorticity; it corresponds to the smaller spacing 
in Z. 

After t — 15, the resolution demands increase 
rapidly. The results at t = 18.75 and 22.50 from 
a properly resolved calculation are displayed in fig- 
ure 10. The computed flow field displays no grid- 
scale oscillations at these times. These results were 
obtained from a calculation on a 64 3 grid until t = 15; 
then the grid was refined to 96 3 , refined further to 
96 x 128 X 128 at t = 16.88, refined yet again to 
96 x 128 x 162 at t = 18.28, refined to 96 x 128 x 192 
at t = 19.22, and refined to 96x 128x216 at t = 21.10. 
The grid was refined by Fourier-Chebyshev interpo- 
lation of the primitive variables from the coarser grid 
onto the finer one. 
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The Fourier-Chebyshev coefficients of the solution 

reveal when additional resolution is required. We 

used the quantities E^ x {t) = ma xEk x ^ z (t), E n (t), 

k z 

and Ek z (t) = ma xE^ x ^ z {t) to assess the resolution 
k x ’ 

requirements in the X-, F-, and Z-directions, respec- 
tively. Several such energy spectra are furnished in 
figure 11. The circles labeled X, for example, are 
a plot of Efr x as a function of the streamwise mode 
k x . Of these examples, only the 64 3 grid results at 
t = 18.75 represent an inadequately resolved calcu- 
lation. The tail of the spanwise spectrum is only 
four orders of magnitude lower than the fundamen- 
tal mode. This corresponds to a truncation error 
in the velocities themselves on the order of 1 per- 
cent. The consequences are clearly apparent in fig- 
ure 9. In this instance, the resolution in F is also 
inadequate, with truncation errors on the order of 
0.1 percent. In no other case shown in figure 11 does 
the truncation error exceed 0.01 percent, that is, an 
eight-order-of-magnitude drop in the energy. These 
computations were performed in 32-bit arithmetic. 
Hence, one should expect no better than five-digit 
accuracy. 

The grid refinements cited above were performed 
when the tail of the energy spectrum in any direction 
reached the 10“ 11 to 10“ 12 level. Moreover, refine- 
ments were performed in only those directions with 
spectral tails this large. For this transition simula- 
tion, the spanwise direction is clearly the most de- 
manding because of the narrow spanwise dimensions 
of the vortex loop near its tip. 

The 32 x 50 x 32 FFDM results for t = 22.5 
are presented in figure 12. Very little of the struc- 
ture of the detached shear layer has been retained. 
The streamwise vorticity plot indicates why this 
computed flow field deserves to be characterized as 
“numerical turbulence” . 

Figure 13 gives the time evolution of several 
Fourier harmonics in the two FFDM calculations. 
The coarse grid fails to capture the full strength of 
the secondary instability, the principal component of 
which is the (1,1) mode. This is consistent with the 
earlier result for linear modes that the coarse grid 
calculation damped these modes excessively. Note 
that the energy contained in the (1,1) mode rises no- 
ticeably less rapidly in the coarse grid simulations. 

The figures given in this section make it appar- 
ent that the use of a coarse finite-difference method 
in the vertical direction instead of either a fine-grid 
finite-difference method or a Chebyshev method pro- 
duces results that differ visibly from the correct solu- 
tion. Moreover, inadequate resolution in the Fourier 
directions clearly leads to numerical turbulence. 


The 32x50x32 FFDM calculation was performed 
to furnish a direct comparison of the present, highly 
resolved simulation with the results of 
Biringen (1984). He continued his simulation on this 
coarse grid all the way to ^ = 80. The results pre- 
sented here indicate that nonphysical, numerical ef- 
fects already invalidate such a calculation by t = 20. 

We close this section with some comments on 
aliasing. Most reported transition calculations have 
been de-aliased, at least in the horizontal directions. 
(The only exceptions appear to be the results of 
Marcus (1984) and of Zang and Hussaini (1985a, 
1985b).) Yet those rigorous numerical analysis re- 
sults which are presently available all suggest that 
de-aliasing is not necessary. (Of course, the differen- 
tial equations must be used in a form which guaran- 
tees conservation of some quadratic integral in ideal 
flow.) For a time-dependent, one-dimensional, lin- 
ear model problem, Kriess and Oliger (1979) have 
shown that the aliasing error is of the same order as 
the truncation error. More importantly, Maday and 
Quarteroni (1982) have proven that aliased and 
de-aliased Fourier spectral approximations to the 
steady, three-dimensional, fully periodic Navier- 
Stokes equations have the same asymptotic rate of 
convergence. Our conclusion is that if aliasing is 
a problem, then the resolution of the calculation is 
marginal anyhow. 

This conclusion is supported by the results given 
in figure 14 for a horizontally de-aliased version of the 
FCBM calculation on a 48 x 64 x 48 grid. The highest 
16 modes in each direction were forced to zero after 
every time step. This removes horizontal aliasing 
interactions (Orszag 1971a). The results here at 
t = 15 are slightly better than the corresponding 
results in figure 8. The de-aliased results at t = 18.75 
are not reliable, but this is a problem of resolution. 
Figure 15 compares the resolution at t = 15. 

All our experience with respect to de-aliasing of 
transition calculations may be summarized as fol- 
lows: Both aliased and de-aliased calculations are 
valid until they lose resolution; the aliased calcula- 
tion loses resolution slightly sooner than a de-aliased 
calculation with an equal number of active modes. 

This experience is contrary to several opinions 
which are cited in the literature. The work of 
Schamel and Elsasser (1976) is frequently invoked 
by those averse to aliased calculations. We have 
repeated one of their examples (the Korteweg-de 
Vries equation) and do not reproduce their results. 
In fact, we find virtually no difference between the 
aliased and de-aliased results. Moreover, Fornberg 
and Whitham (1978) experienced no difficulty with 
aliased calculations of the Korteweg-de Vries equa- 
tion. Moser, Moin, and Leonard (1983) claimed that 
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aliased calculations of Taylor-Couette flow give unac- 
ceptable results. However, the only conclusion that is 
warranted by their numerical examples is that a well- 
resolved de-aliased calculation is superior to a poorly 
resolved aliased calculation. (Their poorly resolved 
de-aliased calculation is just as bad as their sin- 
gle, poorly resolved aliased calculation.) Of course, 
Marcus (1984) performed aliased calculations for this 
problem with no discernible difficulty. 

Hairpin Vortex Structure 

This section is devoted to the detailed structure 
of the hairpin vortex. This is the principal (un- 
steady) intermediate state between the emergence of 
the lambda vortices and the first turbulent activity 
in the transition process. Numerical simulations have 
been performed at the supercritical Reynolds number 
of 8000 as well as at the subcritical Reynolds number 
of 1500. Both K- and H-type transitions have been 
computed. In all cases, the fundamental wave num- 
bers a and j3 were unity. Table II lists the temporal 
eigenfrequencies of the least stable two-dimensional 
TS wave, three-dimensional TS wave, and subhar- 
monic Squire mode under these conditions. 

For each calculation, the grid was refined as 
needed to maintain a truncation error of less than 
0.01 percent in the primitive variables. Some addi- 
tional information on the vertical resolution is pro- 
vided in table III. It lists the location of the (bot- 
tom) critical layer of the two-dimensional TS wave 
and the edge of viscous sublayer as estimated by 
Y v = — 1 + where Y w = vfu T — \j\flR for the 
unpertu rbed lamin ar flow. The wall friction velocity 
u T is \Jv duo/dY. The number of points between 
the wall and the critical layer is N c ; hence, N c /N y 
denotes the fraction of the vertical grid in this region. 
For the viscous sublayer, N u is defined similarly. 

K-Type Vortex at R = 1500 

Subcritical K-type transition at R — 1500 was 
discussed at length in the preceding section. The 
initial TS wave amplitudes were ~ 0.11 and 
e SD ~ 0.10. These initial amplitudes, chosen to 
match those of Biringen (1984), are certainly quite 
large. However, at this low Reynolds number they 
need to be this large in order to trigger the formation 
of a hairpin vortex. 

The time history of the kinetic energy in the low- 
est Fourier harmonics is indicated in figure 16. In the 
early stages (t < 5), the two- and three-dimensional 
TS waves decay linearly. Meanwhile, the flow field 
develops appreciable higher order components which 
are generated by nonlinear effects. The secondary 
instability is evident in the subsequent growth of the 


(1,1) mode. By t = 18.75, the flow has achieved the 
state illustrated in figure 4. This highly developed 
hairpin vortex (the two upper surfaces in fig. 4(a)) 
has evolved far past the stage which can be modeled 
by Floquet theory. 

Although three-dimensional plots such as figure 4 
readily depict the general flow field features, they 
are not capable of conveying detailed information. 
For this, we must resort to a collection of two- 
dimensional cross sections. The vertical shear and 
the perturbation velocities in the peak plane are dis- 
played in figures 17 and 18. The velocity vectors are 
shown, for clarity, at only a subset of the computa- 
tional grid. The full channel is displayed in figure 18 
to illustrate the degree to which the shear layers in 
the two halves of the channel are interacting. By 
these late times, distinct kinks have developed in the 
shear layer. Some streamwise velocity profiles are 
given in figure 19. The detached shear layers clearly 
correspond to inflectional velocity profiles. Williams, 
Fasel, and Hama (1984) have argued that this shear 
is created as the fluid surrounding the hairpin vortex 
moves up and over it. 

The evolution of the vertical shear is similar until 
t = 22.5 to that reported by Kovasznay, Komoda, 
and Vasudeva (1962) for controlled transition in a 
boundary layer. At about t = 22.5 (fig. 18), the shear 
layers which originate at the bottom and the top of 
the channel begin to overlap. This process does not 
occur in controlled transition in the boundary layer 
or in the channel at higher Reynolds numbers. 

Figures 20 to 23 present detailed information 
on the hairpin vortex at t — 18.75. Its shape is 
most apparent in the plots of streamwise vorticity 
(fig. 20(a)). Solid lines denote fluid rotating clock- 
wise, and dashed lines denote fluid rotating counter- 
clockwise. The legs of the vortex appear to emerge 
near X/L x = 1/8 in the vicinity of the critical layer 
of the two-dimensional TS wave. Since the spanwise 
direction has been compressed by a factor of 4.19 in 
these plots, the legs of the vortex are even flatter 
than they appear in figure 20(a). As X increases, 
the vortex becomes less flat, lifts away from the wall, 
and pinches toward the peak plane. 

Two regions of enhanced spanwise vorticity 
(fig. 20(b)) can be identified: the detached shear 
layer, which is just above the tip of the hairpin vor- 
tex, and the region between the legs of the vortex and 
the wall. There is also a local minimum in the (abso- 
lute value of the) spanwise vorticity just underneath 
the detached shear layer. Thus, the hairpin vortex is 
characterized by a region of large streamwise vortic- 
ity and small spanwise vorticity. These are the same 
general features which were measured by Williams, 
Fasel, and Hama (1984) and observed numerically by 
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Zang and Hussaini (1985b) for the controlled bound- 
ary layer. 

Figure 21, which displays velocity vectors in Y- Z 
cross sections, reinforces the picture of the hairpin 
vortex suggested by figure 20. Figure 22 does the 
same. It presents contour plots of the helicity, defined 

by 

h = u * u (9) 

Helicity has been proposed as one indication of co- 
herent structures in turbulent flows (Levich and 
Tsinober 1983). (The universality of this has, 
however, been criticized by Speziale (1987) on the 
grounds that helicity is not a Galilean invariant quan- 
tity.) This quantity provides a less cluttered picture 
of the streamwise vorticity. 

The effects of vortex stretching and diffusion are 
illustrated in figure 23. The equations of motion 
imply that 

cS J \ ^ — J ^ ^ ^ + J : Vu; dV 

< 10 ) 

The integrand of the first term on the right side mea- 
sures the rate of vortex stretching and the second 
term measures the rate of vorticity diffusion. The 
surface integral accounts for the generation of vortic- 
ity at the walls. The first two integrands are plotted 
in figures 23(a) and 23(b). The most rapid stretching 
is clearly associated with the regions of most intense 
spanwise vorticity, while the diffusion is most signif- 
icant at the edges of these regions. The net change 
in vorticity is dominated by stretching effects. 

Figures 24 to 26 supply the vorticity contours at 
t = 15, 22.5, and 25.0. Together with figure 20, 
they provide a time history of the hairpin vortex 
and the detached shear layer. At all four of these 
times, the shear layer does indeed lie immediately 
above the vortex tip. Hence, it may be used to 
estimate the speed at which the vortex is traveling 
in the streamwise direction. Table IV summarizes 
information extracted from figures 17 and 18. The 
location of the shear layer in the peak plane is 
denoted by (Xsl,Tsl). They are the coordinates 
of the center of the shear layer. The shear layer 
speed wsl i s calculated from this information, after 
taking into account the speed at which the coordinate 
system itself is moving (0.36 in this case). The 
speed of the unperturbed mean flow at the height 
of the shear layer is given for reference. The vortex 
moves faster as it rises away from the wall and slows 
down after it crosses the center of the channel. It is 
indeed moving slower than the surrounding fluid, as 
suggested by Williams, Fasel, and Hama (1984). 


By t = 22.5, (fig. 25) a counterrotating vortex has 
emerged inside the hairpin vortex. At X/L x = 1/2 
and 3/4 the spanwise vorticity inside the neck of the 
hairpin vortex is so small that the streamwise vor- 
ticity is dominant. Hence, this is a second region, 
in addition to the hairpin vortex itself, where the 
streamwise vorticity is large and the spanwise vortic- 
ity is small. This feature is still present at t = 25. 
At this stage, it appears that the neck of the hairpin 
vortex is in the process of splitting into two parts. 
The appendix contains additional data on the flow 
field at these times. 

H-Type Vortex at R = 1500 

The initial conditions for the H-type transition 
at R = 1500 again were t2D — 0.11 and €%£> = 0.10, 
but the three-dimensional wave was the subharmonic 
Squire mode (with s x = 2). The simulation started 
on a 64 x 64 x 32 grid and was refined as needed — by 
t = 32.5, a 96 x 96 x 162 grid was in use. 

The energy history is provided in figure 27. The 
initial three-dimensional wave is denoted by (1/2,1), 
indicating that its streamwise wave number is only 
half that of the two-dimensional TS wave. Both 
modes decay linearly in the early stages, but by t = 8 
the secondary instability has asserted itself and is 
reflected by the growth of the (1/2,1) mode. 

The structure at t = 30 is typical of the later 
stages. Figure 28 gives a three-dimensional picture 
of the vorticity. The vortices are, of course, arranged 
in the familiar staggered pattern. The vertical shear 
in a peak (and one other) plane is given in figure 29. 
This shear layer is roughly equal in strength and 
height to the K-type one at t = 15. (Note that the 
X-axis has been compressed by a factor of 2 relative 
to the K-type plots.) Vorticity contours are provided 
in figure 30. The structures here are similar to those 
for the K-type vortex, albeit more complicated due 
to the substantial interference effects of neighboring 
vortices. Helicity, stretching, and diffusion contours 
are supplied in the appendix. 

The hairpin vortex evolves much less rapidly for 
the subharmonic instability than for the fundamental 
one. The vorticity structure is qualitatively similar 
and only slightly more intense at t = 40 than at 
t = 30. The K-type vortex, however, undergoes a 
very rapid evolution between t = 15 and t = 25. 
Indeed, the relative weakness of the subharmonic 
vortex system is responsible for its less stringent 
resolution requirements. 

K-Type Vortex at R = 8000 

The next two simulations were performed under 
conditions much closer to those of actual transition — 
the Reynolds number was supercritical ( R ~ 8000), 
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the initial two-dimensional TS wave began at the 
2-percent level and the initial three-dimensional wave 
had the very low amplitude of 0.02 percent. The res- 
olution requirements, especially in the normal direc- 
tion, were more severe at the higher Reynolds num- 
ber. Hence, the simulations were carried to a less 
advanced state than those for R = 1500. 

The harmonic energy history of the K-type tran- 
sition is furnished in figure 31. This simulation was 
conducted for more than seven periods of the two- 
dimensional TS wave. For at least the first six pe- 
riods (t < 153), the growth of the two-dimensional 
wave is described very well by linear theory. The 
secondary instability induces the growth of the three- 
dimensional wave after one period, but because of the 
low initial amplitude of the three-dimensional dis- 
turbance, the hairpin vortex does not emerge un- 
til roughly 6 1/2 periods ( t > 165). The subse- 
quent evolution is extremely rapid. The simulation 
was halted at t — 182.21 (7 1/6 periods) when the 
72 x 162 x 192 grid had just become inadequate. 

Several illustrations of the resulting structure at 
t = 182.21 are provided here in figures 32 and 33. A 
comparison of figure 32 with figure 17 and of figure 33 
with figures 20, 24, 25, and 26 suggests that the 
t = 18.75 state of the low Reynolds number vortex 
is the closest match to the R = 8000 vortex. (In 
comparing these figures, note that the R = 8000 
vortex is shifted forward by 1/8 L x in the streamwise 
direction.) 

Clearly, the only substantial difference between 
these two vortices are in scale and intensity. The 
R = 8000 structures are narrower and closer to 
the wall. Despite streamwise velocity perturbations 
which are 20 percent smaller, the detached shear 
layer is 100 percent more intense for R = 8000 than 
for R = 1500. This is a result of larger gradients at 
the higher Reynolds number. The streamwise vor- 
ticities for the two cases, however, have comparable 
amplitudes. Helicity and vortex stretching plots are 
supplied in the appendix. 

H-Type Vortex at R = 8000 

The harmonic energy history of the H-type tran- 
sition at R = 8000 is supplied in figure 34. Like its 
R — 1500 counterpart, the instability is slightly 
weaker than the K-type secondary instability, judged 
by the rate at which grid refinements were re- 
quired. This simulation was halted at t = 186.45 
(7 1/3 periods). The final grid was 64 x 128 x 128. 

The structure of this vortex is illustrated in fig- 
ures 35 and 36. It resembles its lower Reynolds num- 
ber counterpart in all essential respects but one: the 
vortices are less intertwined at R = 8000. 


Discussion 

The structure of the hairpin vortex is difficult 
to resolve, even at low Reynolds numbers. A de- 
tailed description of the resolution requirements of 
one case was reported at length in an earlier section. 
These results, and similar studies for the other three 
cases, lead us to the following conclusions. (1) The 
streamwise direction requires by far the least resolu- 
tion. Only in the later stages of transition, when the 
detached shear layer develops kinks, does high reso- 
lution in this direction become necessary. (2) Even 
at low Reynolds number, the spanwise direction is 
very demanding. (3) The resolution in the vertical 
direction is the most sensitive to Reynolds number. 
(4) A grid refinement is needed in any direction when 
the tail of the energy spectrum reaches 10~ 8 of the 
low-frequency value. 

These remarks are buttressed by figure 37, which 
presents resolution plots for each of the four simu- 
lations at the last stages reported here. The only 
case with substantial high-frequency content in the 
X-direction is the R = 1500, K-type vortex at 
t = 25. This is the only simulation in which there 
are at least two streamwise kinks. At the slightly ear- 
lier time of t — 22.5, the second kink is absent and 
so is the streamwise high-frequency content (fig. 11). 
The high spanwise resolution is needed because of the 
sharp gradient and complex structures which develop 
at the tip of the vortex. 

The difficulties presented by higher Reynolds 
number calculations are exemplified by the 
R = 8000, K-type vortex. The vertical frequency 
spectrum (fig. 37) has a very long tail at t = 182.21. 
In fact, this direction fails the eight-decade guide- 
line. Indeed, this marginal lack of resolution is 
responsible for the oscillations in the 0.2 vertical 
shear contour underneath the detached shear layer 
in figure 32. These calculations were all performed 
on standard Chebyshev grids (Yy = — cos ( 7 xj/N y ), 
j = 0, 1, . . . , N y ). Numerical simulations such as 
these would surely benefit from a more highly tuned 
vertical grid. A mapping can be readily incorporated 
into both spectral algorithms used for these calcu- 
lations (Zang and Hussaini 1985a and 1986) for at 
most a 50-percent increase in running time. Map- 
pings do, however, make algorithms which rely on 
direct solution techniques for special basis functions 
(Moser, Moin, and Leonard 1983) impractical. More- 
over, they necessitate the more expensive matrix di- 
agonalization versions of the Orszag-Kells (1980) and 
Kleiser-Schumann (1984) methods. 

Another useful modification to the algorithm 
would involve using a nonuniform grid in the 
Z-direction that clustered the points near the peak 
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plane. This would complicate the solution of the im- 
plicit equations for pressure and for normal diffusion. 
The spectral multigrid approach suggested by Zang 
and Hussaini (1986) would be able to handle this 
complication. 

The K- and H-type transitions produce similar 
hairpin vortices, although the K-type vortices appear 
to be more intense. Conceivably, the H-type vortices 
are stronger at different spanwise wave numbers. 

The vortices which arise at higher Reynolds num- 
bers are less diffuse and closer to the wall than those 
at lower Reynolds numbers. Vortex stretching ef- 
fects predominate over vortex diffusion in most of the 
hairpin vortex. The streamwise and spanwise vortic- 
ity patterns conform to the large-scale measurements 
of Williams, Fasel, and Hama (1984) for a K-type 
boundary-layer transition. Examination of the vor- 
ticity patterns displayed in this report reveals that 
the detached shear layer is in fact two separate shear 
layers above the legs of the vortex which merge into 
a single shear layer above the vortex tip. 


The most surprising result was the robustness 
of the time-dependent hairpin vortex. None of the 
transitional flows simulated here became turbulent 
before all the available resolution (nearly 3 million 
grid points) was exhausted. There were no signs of 
violent inflectional instabilities emanating from the 
detached shear layer. (Several additional runs of 
the R = 1500, K-type vortex included a random 
component in the initial velocity field. These were 
damped very rapidly by the high viscosity of this 
flow. Even for an initial streamwise root-mean- 
square amplitude of 0.2 percent, these fluctuations 
decayed so much that by t = 18.75 they could no 
longer be detected in flow-field plots. Evidently 
larger amplitudes and higher Reynolds numbers are 
required for these instabilities to occur.) Seemingly 
chaotic behavior was achievable only by resorting to 
an obviously inadequate grid. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
December 18, 1986 
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Appendix 

Supplementary Illustrations 

This appendix is directed toward those interested 
in additional details from these four simulations. 

Figures 38 to 43 pertain to the R = 1500, 
K-type transition. The first three present a more 
complete picture of the vorticity fields at t = 15, 
22.5, and 25. The last case is especially complex. 
Figure 41 contains cross-stream velocity vectors at 
t = 22.5. The last four frames indicate that the 
hairpin vortex actually dips back down to the wall 
(near X/L x = 5/8) just before the tip. Figure 42 
presents the helicity at t = 22.5. Figure 43 supplies 


not only the vortex stretching and diffusion rates at 
t = 22.5 but also their sum. (Because of the intensity 
of the stretching, the contours in the most intense 
regions have been omitted.) The net change of vor- 
ticity is conveyed by the last part of figure 43. Note 
that between X/L x = 0 and X/L x = 1/2, there is a 
substantial region along the wall in which there is a 
net decrease of vorticity. Elsewhere, however, vortex 
stretching dominates the diffusion effects. 

Figures 44 to 49 are self-explanatory and pertain 
to the other three simulations. Vortex diffusion plots 
have not been supplied for the R = 8000 cases 
because of plotting difficulties arising from very large 
gradients near the wall. 
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Nomenclature 

a vertical grid stretching parameter 

Ek x ,k z kinetic energy in (k X) k z ) Fourier 
component 

E n kinetic energy in nth Chebyshev 

component 

h dimensional channel half-width 

k x scaled streamwise Fourier harmonic 

k z scaled spanwise Fourier harmonic 

k x streamwise Fourier harmonic 

k z spanwise Fourier harmonic 

L x periodicity length in X-direction 

L z periodicity length in Z-direction 

N c number of grid points inside critical 

layer 

N x number of grid points in X-direction 

N y number of grid points in Y-direction 

N z number of grid points in Z-direction 

N u number of grid points inside viscous 

sublayer 

P total pressure 

p pressure 

R Reynolds number, u c h/v 

s x number of subharmonics in 

X-direction 

s z number of subharmonics in 

Z-direction 

T n Chebyshev polynomial of degree n 

t time 

u total velocity vector 

u kxjkz velocity after transforms in X and Z 
u n velocity after transform in Y 

u streamwise velocity component 

u c dimensional centerline velocity 

^SL mean flow at shear layer height 


uq mean flow 

U2 d two-dimensional eigenfunction 

U3£> three-dimensional eigenfunction 

v normal velocity component 

w spanwise velocity component 

X streamwise coordinate 

Xsl X position of shear layer 

Y normal coordinate 

Y c Y position of critical layer of two- 

dimensional TS wave 

Y sl Y position of shear layer 

Y u Y position of edge of viscous sublayer 

Z spanwise coordinate 

a streamwise wave number 

/? spanwise wave number 

6 ^/ Kronecker delta function 

A Z grid spacing in spanwise direction 

62 d initial amplitude of two-dimensional 

TS wave 

63/) initial amplitude of three-dimensional 

TS wave 

v dimensional kinematic viscosity 

£ computational coordinate in normal 

direction 

lj vorticity 

o^sqr temporal frequency of subharmonic 

squire mode 

0J2D temporal frequency of two-dimensional 

TS wave 

u?3£> temporal frequency of three- 

dimensional TS wave 

Abbreviations: 

FCBM Fourier-Chebyshev method 

FFDM Fourier finite-difference method 

TS Tollmien-Schlichting 
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Table I. Effect of Vertical Resolution on Growth Rate 
[R = 1500, three-dimensional TS wave] 


Algorithm 

Ny 

Calculated 

Im(w 3 £>) 

Error in 
Im(w 3 £>) 

Exact 


-0.028230 


FCBM 

64 

-.028226 

-0.000004 

FFDM 

50 

-.032030 

-.003800 

FFDM 

100 

-.029230 

-.001000 

FFDM 

200 

-.028482 

-.000252 


Table II. Eigenfrequencies of the Least Stable Linear Waves 


R 

w 2D 

W 3D 

w Sqr 

1500 

8000 

0.326299 - 0.028206i 
0.247075 + 0.002664i 

0.401293 - 0.028230i 
0.292106 - 0.022388i 

0.171386 — 0.0926171 
0.099426 - 0.0553271 


Table III. Resolution Information 


R 


N e /Ny 

Y u 

N u /Ny 

1500 

-0.821 

0.193 

-0.982 

0.061 

8000 

-.868 

.166 

-.992 

.040 


Table IV. Position and Speed of Detached Shear Layer for the R — 1500, K-Type Transition 


t 

X sl/Lx 

^SL 

U SL 

®o(^sl) 

15.000 

0.56 

-0.31 

0.69 

0.90 

18.750 

.78 

-.16 

.78 

.98 

20.625 

.92 

-.06 

.83 

1.00 

22.500 

1.06 

.12 

.81 

.98 

25.000 

1.23 

.31 

.79 

.90 
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Figure 1. Lambda vortex patterns in transitional flows. 



Figure 2. Channel flow geometry. The streamwise coordinate is X, the spanwise Z, and the normal (or vertica 
coordinate is Y. 
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2-D TS WAVE 

(k x ,k z ) = (1,0) 



3-D TS WAVE 

(k x ,k z ) = (1,-1) 



SPANWISE MODE 

(k x ,k z ) = (0,1) 



Figure 3. Horizontal structure of several Fourier components. The flow direction is indicated by the arrows. 
Two fundamental wavelengths are plotted in both streamwise and spanwise directions. 
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(a) Streamwise vorticity component. 



L-86-412 


(b) Spanwise vorticity component. 

Figure 4. Vorticity components of the R = 1500, K-type transition at t = 18.75. 
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32 x 200 x 32 



Figure 5. Vertical shear ( du/dY ) in the peak plane ( Z/L z — 1) at t = 7.5 from two FFDM calculations. 



Figure 6. Vertical shear in the peak plane (left) and streamwise vorticity at X/L x — 1/2 (right) from a FFDM 
calculation on a 32 x 50 x 32 grid at t — 15. 






Figure 7. Vertical shear in the peak plane (left) and streamwise vorticity at X/L x = 1/2 from a FFDM 
calculation on a 32 x 200 x 32 grid at t = 15. 
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Figure 8. Vertical shear in the peak plane (left) and streamwise vorticity near the tip of the vortex loop (right) 
from a FCBM calculation on a 32 X 64 X 32 grid. 















Mode Mode 

Figure 11. Energy spectra in each of the coordinate directions for the FCBM calculations. 
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Mode Mode 

Figure 15. Energy spectra in each of the coordinate directions for a 32 x 64 x 32 FCBM calculation and 
horizontally de-aliased 48 x 64 x 48 FCBM calculation. 
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Figure 16. Kinetic energy in several (k x ,k z ) harmonics for the R = 1500, K-type transition. 
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Figure 21. Cross-stream (v-w) velocity vectors for the R = 1500, K-type transition at t = 18.75. 
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(b) Vortex diffusion rate 
Figure 23. Concluded. 
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(a) Streamwise vorticity component. 



L-86-414 
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(b) Spanwise vorticity component. 

Figure 28. Vorticity components of the R = 1500, H-type transition at t = 30. 































(a) Streamwise vorticity. 
Figure 33. Vorticity contours for the R = 8000, K-ty 

















(b) Spanwise vorticity. 
Figure 33. Concluded. 
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Time 

Figure 34. Kinetic energy in several (k x ,k z ) harmonics for the R — 8000, H-type transition. 
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Figure 35. Vertical shear at Z/L z — 1/4 and in one peak plane ( Z/L z = 1/2) for the R — 8000, H-type 
transition at t = 184.45. 







(a) Streamwise vorticity. 

Figure 36. Vorticity contours for the R = 8000, H-type transition at t — 184.45. 










(b) Spanwise vorticity. 
Figure 36. Concluded. 















Figure 37. Energy spectra for 




Mode 


four transition simulations. 























(a) Stre; 


Figure 39. Vorticity contours for 














(b) Spanwise vorticity. 
Figure 39. Concluded. 
















(a) Streamwise vorticity. 

Figure 40. Vorticity contours for the R = 1500, K-type transition at t 














(b) Spanwise vorticity. 
Figure 40. Concluded. 
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Figure 41. Cross-stream (v-w) velocity vectors for the R = 1500, K-type transition at t = 22.5. 
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